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1.   Introduction 

In  recent  years  there  has  appeared  in  the  literature  a 
surge  in  the  number  of  papers  dealing  with  numerical  solutions 
of  partial  differential  equations.   And,  usually,  the  difference 
methods  employed  are  of  first  or  second  order  of  accuracy. 
This  restriction  is  not  an  arbitrary  one,  but  rather,  is  related 
to  the  fact  that  computing  machines  have  been  relatively  slow 
and  their  high  speed  memory  capacity  has  been  small;  hence  a 
usable  computational  scheme  must  necessarily  have  the  attribute 
of  simplicity.   In  problems  of  more  than  one  space  dimension, 
even  greater  emphasis  is  placed  on  simplicity. 

It  is  anticipated,  however,  that  a  new  era  of  computability 
is  almost  upon  us.   We  are  referring  to  the  use  of  parallel 
processors,  i.e.  N-serial  type  computing  processors,  each  of 
which  is   synchronized  and  each  of  which  can  communicate  with 
the  other  processors  through  a  common  memory  or  central 

r  o 

controller.   The  value  of  N  may  be   from    2   to  2   and  the 
arithmetic  speed  of  each  individual  processing  unit  will  be 
in  the  sub-microsecond  range.   By  proper  organizing  of  the  data, 
each  mesh  point  or  string  of  mesh  points   may  have  its  own 
central  processor,  which  means  the  solution  on  the  entire  mesh 
may  be  advanced  essentially  simultaneously.   For  such  a  class 
of  computing  machines  the  requirement  of  simplicity  for  the 
difference  scheme  may  be  relaxed. 

In  this  note,  we  propose  a  class  of  difference  schemes  for 
hyperbolic  problems  in  one  and  two  space  dimensions.   The 
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methods  are  applicable  to  nonlinear  initial  value  problems; 

they  are  uniformly  third  order  accurate  on  both  the  space  variables 

and  time  and  are  similar  to  methods  proposed  by  Strang  [8]. 

2 .   A  Third  Order  Difference  Operator 

We  will  construct  a  difference  operator  which  is  uniformly 
accurate  to  third  order  in  each  of  the  space  and  time  increments 
Ax  and  At.   There  are  schemes  which  are  third  or  fourth  order 
accurate  in  the  space  variable,  but  second  order  in  the  time  step. 
Such  methods  are,  therefore,  not  uniformly  accurate  in  both  the 
independent  variables.   The  approximation  scheme  described  here 
is  constructed  in  divergent  form  --  Just  as  the  original  differ- 
ential equation  is  written  in  divergent  or  conservation  form. 
In  order  to  describe  the  derivation  consider  the  differential 
equation  in  one  dimension 

u   =  F(u,x,t,u  ) 
(2.1)  ^  ^ 

u(x,0)    =    U„(x)        _oo  <  X  <  °° 

in  which  the  flux  is  computed  by  evaluating  F.   Of  course,  for 
partial  differential  equations,  F  cannot  be  computed  exactly 
since  it  depends  on  derivatives  in  the  space  variables.   An 
approximate  evaluation  of  F  can  be  obtained  if  the  space 
derivative,  u  ,  is  replaced  by  a  difference  approximation  6u. 
In  this  paper  we  look  at  Runge-Kutta  type  approximations  to 
(2.1)  which  are  third  order  accurate  and  for  which  the  algorithm 

u(t+At)  =  u(t)  +  ^  (k^+^k2+ko) 
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where 

k-^^  =  At  F(t^,u^) 

k2  =  At  F(t.+At/2,  u^+k^/2) 

k^  =  At  F(t^+At  ,   u^-k-^+2k2) 

Is  a  third  order  Runge-Kutta  method  of  integration  for  first 
order  ordinary  differential  equations  and  is  analogous  to  the 
method  derived  for  partial  differential  equations.   Since  three 
evaluations  of  F(x,t)  are  required,  we  would  expect  the  same 
number  of  evaluations  of  difference  approximations  to  the  flux 
for  the  partial  differential  equation  to  be  required.   The 
form  of  the  function  F  for  the  partial  differential  equation 

(2.2)  u^  =  f^ 

leads  to  the  requirement  that  the  approximation  to  F,  F  differ 
only  by  terms  of  0(At  )  so  that  if  u  is  a  difference  approxima- 
tion satisfying 

u(t+At)  =  u(t)  +  At  P 

then  |u(t)-u(t)|  =  0(At  ).  Here 

•5 

F(x.,t  ,u.,6u,...)  =  F(x,t ,u,u  ,  .  .  .  )  +  0(At  )  , 
1   n   1  x 

The  great  advantage  of  Runge-Kutta  methods  (one  of  which  is 
Richtmyer's  Lax-Wendroff  two  step  method  [1])  is  that  to  achieve 
higher  order  accuracy  in  approximating  u(t)  only  repeated 
evaluations  of  F  are  required.   For  complicated  F,  however. 
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this  advantage  becomes  a  disadvantage  so  that  it  may  be  desirable 
to  use  methods  which  do  not  require  multiple  evaluations  of  F(x) 
having  only  slight  changes  in  the  value  of  the  argument,  x,  of  F. 

As  In  ordinary  differential  equation  theory  it  is  possible 
to  construct  approximation  techniques  to  partial  differential 
equations  using  Taylor  series  methods.   For  high  accuracy  this 
would  require  evaluation  of  higher  derivatives  of  (2,1)  and  for 
coupled  systems  of  equations,  i.e.,  for  systems  of  the  form  (2.2), 
these  derivatives  become  more  and  more  complicated  to  evaluate. 
Lax  and  Wendroff  [2]  showed,  however,  that   for  systems  of  conserva- 
tion laws   given  by  (2.2),  Taylor's  method  can  be  used  to  construct 
an  elegant  second  order  accurate  algorithm.   They  showed  that  u 

X/  X/ 

can  be  evaluated  by  using  the  original  differential  equation  (2.2). 
If  the  matrix  A  =  9f/8u   is  Introduced,  then  (2.2)  can  be  written 
by  using  the  chain  rule  as 


(2.2a) 


u,  =  A(u)u 
t  X 


so  that 


(2.3) 


U^.^.  =  f^ 
tt     tx 


(|^  1^)   =  (Af  ) 
9u  3t  Y       XX 


is  found  in  terms  of  space  derivatives  only.   If  one  wished  to 
construct  a  third  order  method,  it  would  be  necessary  to  compute 
the  time  derivative  of  (2.3)  for  the  next  term  in  the  Taylor 
series.   Unfortunately,  it  is  not  possible  to  then  eliminate  terms 
containing  3A/8t  =  A  so  that  terms  containing  only  space 
derivatives  remain.   It  is  clear  that  the  dependence 
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A  =  A(u,u^) 

will  occur  through  time  differentiation  of  each  term  of  A  and 
back  substitution  for  the  time  derivatives  of  u  from  the  right 
hand  side  of  (2.2)  (where  the  function  P  is  relatively  simple). 
For  the  equations  of  gas  dynamics  this  procedure  will  result 
in  unnecessarily  complicated  difference  methods  with  the 
associated  disadvantage  that  the  form  of  the  algorithm  will 
not  be  conservative.* 

Instead  of  pursuing  this  approach,  we  use  the  alternate 
procedure   which  we  first  discussed,  and  which  was  first  pointed 
out  by  Robert  Richtmyer.   He  showed  that  the  Lax-Wencroff  method 
could  be  written  in  two  steps.   For  each  step,  only  an  evaluation 
of  f   is  required  —  Just  what  one  would  expect  from  a  Runge-Kutta 

A. 

type  method.   The  third  order  method  which  we  now  describe  was 

first  proposed  by  Rusanov,  and  we  repeat  some  of  the  results  that 

are  contained  in  his  paper  [3].   We  consider  a  sequence  of 

iterates  to  the  solution  u(t).   The  r-th  iterate  defines  an 

approximation  to  (2.1),  and  is  given  by 

r-1 
u   =  u   +  At  I      a        F(u   t  )  ,    r  =  1,2,. ..,R. 
s  =  0 

The  function  F  is  evaluated  at  time  t   =  t^  +  t   At  ,  s  =  1 ,2 , . . .  ,R , 

S  0  S        '  5    5        3    5 


IE 

Gideon  Zwas  has  shown,  however,  for  the  scalar  case  where 
A  =  a(u),  that  the  £-th  time  derivative  of  u  can  be  given  by 
the  compact  expression 

u„,  =  (a  u  )    ,  n  =  £-1, 

£t       X  nx  '  ' 

which  preserves  the  conservation  form  of  the  associated 
differential  equation. 
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and  u(0)  =  u^  is  defined  at  t„. 

A  solution  u„  =  u(t-,)  is  obtained  at  t-,  =  t„  +  At.   To 
advance  the  solution  from  t-,  to  tp,  the  above  procedure  is 
repeated  with  t„  replaced  by  t^ .   The  a    are  determined  by 
requiring  that   u„   satisfy  the  Taylor  expansion 

X  A^  f^^^O,    ^    At^  ,^  ^0,  ,  At^  .3  ^0,  _,    ...A, 


up  to  the  required  order  of  accuracy  which  is  three.   The 
quantities  in  parentheses  are  difference  approximations  to  the 
derivatives  of  u^.   To  apply  this  procedure  to  the  partial 

differential  equation  (2.2)  it  is  convenient  to  write  out  the 

(n ) 
sequence  of  iterates  using  the  notation  u.    =  u(x.,t  )  with 
M  &  1        1 '  n 

t   -,  =  t   +  At  and  A  =At/Ax.   We  use  the  following  spatial 
difference  operators  in  the  derivation 


(2.3) 


U  f(x^)  =  I  (^(^1+1/2^  +  ^^^i-1/2^^ 


6  f(x.)  =   f(x.^,/2)  -  f^^-1/2) 


and  thus      y6  f(x.)  =  |  (f(x^_^-|_)  -  f(x._-^)) 

Then  u.    =  u(x.,t   +T,,At)   is  given  by 
1        1 '  n   1         °      "^ 


(2.^)  u.^^/2  =  ^^i+1/2  +  "lO^^^^i+1/2^   ' 

( 2) 
u.    =  u(x. ,  t   +  T„  At)   is  given  by 
1        1 '   n     2         ^  ■^ 

(2.5)       u|2^  -  u|°^  +  a2Q{Xy5fj°h  +  a^^iX^f^^h 
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and  u?"*"   =  u.  '^  =  u(x.,  t  +At)  is  given  by 


n+1        (0)   _^        r,  ,-r  ,  Q     .2s   r^(O)-, 

u.    =  u.    +  aT„{A(I  +  0-,-,6  jyof.  s 


(2.6)  +  a^^{\(l+Q^^&^)    6fj^h 


+  a,^{Ay6f  ^^h 
32   "^   1 


The  sequence  (2.4)  -  (2.6)  is  chosen  to  be  in  flux  divergent 
form  since  the  original  differential  equation  is  in  this  form. 
Equations  (2.4)  and  (2.5)  are  generalizations  of  Richtmyer's 
two  step  method;  it  will  be  shown  that  this  form  will  lead  to 
a  one  parameter  set  of  difference  methods.   Equation  (2.6) 
represents  a  linear  combination  of  central  differences  of  the 
flux   at  the  three  previous  time  levels   0,  x.At  and  TpAt  .  The 

quantity  yu,  ::  ,„  replaces  u.  .,  ,„  for  stability  of  (2.4)  (and 

(2) 
indeed  for    stability  of  u.   ). 

The  prescription  to  find  the  a    is  to  use  operators  (2.3) 

rs 

in  (2.4)-(2.6),  and  then  expand  each  term  in  the  brackets  as  a 

Taylor  series;  for  example 

2 
yf .     =    f .    +   i  d^f      ^^ 


i    ~      i         2      x    i      2 

3 

(2.7)  5f.    =    d    f.    Ax    +    ^  d^f.    ^^ 

1X1  6      X    1      8 


y6f.E    d    f.    Ax   +   J  d^f.    Ax^ 
1        XI  6x1 
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We  use  the  symbol  [e]  to  indicate  that  the  expressions  (2.7) 

4 
are  correct  modulo  terms  (}(Ax  ).   In  order  to  compare  the 

resulting  expression  with  the  third  order  expansion 


2  ^ 

L^^u  ;-^-  +  (d^^^u  )-^ 

(2.8) 


u-^-^l  =    u^  +  (d,u")At  +  (d^.u^^)^  +  (d,,,u^)^  +  0(AtS 


=  u-  -f  (d^f-)At  +  (d.^f'^)^  +  (^ttx^")^  +  ^(^^"^^  ' 


function  f  ^  ^  (t  +T-j^At )     and  f   '^(t„  +  T„At)  must  be  expressed 

as  a  Taylor  expansion  about  f'^  =  f    (t„).   For  instance,  since 

( ?) 
we  want  to  find  y6f.   ,  first  use  (2.7)  and  then  apply  Taylor's 

formula  to  the  result  to  obtain 

2 
y6f{2)  .  (d^f{0))Ax  .  X2At(d^^f^0))Ax  .  ^^^4^  (d^tt^^^^^^ 

+  i    (d3f(0))  Ax3  . 
6x1 


In  a  similar  fashion,  we  obtain  the  required  expansion  for  each 

of  the  bracketed  expressions  in  equation  (2.6).   The  expression 

^     n+1      n  ^-   • 
for  u    results  m 

^n+1  =  ^n  _^  (a^Q  +  a^^^^^x^-*^^  "*■  "32'^2  "^^xt^"*^^^ 

2 
(2.9)  +  ^P  (dxtt^)^t^  -^  ^31^^xx^^^^' 

1  2 

+  -r    (otoA  +  6  ao^e,^  +  a^p)  At  Ax 
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The  terms  containing  a^-,  are  at  most  second  order  accurate  so 

to  preserve  the  accuracy  of  u    ,  a^-,  =  0. 

n+1 
Comparing  (2.8)  with  (2.9)  we  see  that  for  u     to  be  third 

order  accurate,   equation  (2.6)  must  have  coefficients  a^   and 

6^-,  which  satisfy 

Ctnp.   +    d-^^         —    J-  y  Ct  ^  p  T  p    "~   "p"    3 

(2.10)  2 

32  2    1 
0I3Q  +  6  ci^Qe^-^  +  a^p  =  0   ,    — 2 ^  ^  • 

These  relations  Imply  that  a^p  =  3/^,   ot^Q  =  1/^,   Xp  =  2/3 
and  9^-,  =  -2/3.   Again,  by  using  this  expansion  procedure  on 
the  bracketed  expressions  In  equation  (2.5),  we  obtain 

(2.11)  u*^^^  =  u'^  +  a2Q(d^f)  At  +  ap^((d^f)At  +  (d^^f)T^At^)  . 

(2) 
For  the  above  equation  for  u     to  differ  from  the  Taylor 

expansion  for  u     about  u   by  only  terms  of  (}(At  ),  the  a^ 

must  satisfy 

1   2 

(2.12)  a^Q  +  a2i  ~  "^2  '      '^2l'''l  ^  2"  "^2  ' 

Similarly,  for  equation  (2.^)  to  yield  first  order  accurate  data 
for  u^-^\ 

(2.13)  a-j^Q  =  T-^    . 

Hence,  we  have  specified  the  coefficients  to  within  one  parameter, 
namely  t-,  . 
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For  T  -    1/3,  we  have  ap„  =  0  and  a^-.    =    2/3.   In  this  case 
the  difference  equations  become 


(2.1^a) 


u 


(1) 
i+1/2 


1/n   .nN,lr,.„n    ^n^-. 
2  (•^i+l+^i^  +  3  {Mf.+i-f.)] 


(2.14b) 


(2) 
u. 

1 


^i  +  3  ^^^^i+1/2  -  ^i-1/2^^ 


(2.l4c)     u 


n+1 


u 


n  ,  1  rA 


1  ^  i  ^^^-^'1.2'  ^'1,1-  7f?_i-^  2f;_2) 


,  3  A(A2)         f,(2)s. 


These  equations  are  the  analog  of  the  integral  of  the 
conservation  law.   To  see  this  we  use  the  figure  to  define 
integral  quantities : 


U-,  = 


u(x  jt-.  )  dx 


f  (x^,t  jU^) 


\ 


'///////////////m 


f  (x„,t ,u^) 


X^-   T^ 


Uq   = 


uCXjtg)    dx 


Figure    1 


-10- 


Integrate  (2.2)  over  the  shaded  region  In  the  space  time 
domain  to  obtain 


U^  =  Uq  +  J   f(x2,t,U2)  dt  - 


1 
f 


f  (x-,  ,t  ,u-,  )  dt 


In  the  above  expression  u-,(t)  and  ijp(t)  are  the  values  of  u  at 
the  net  points  x,  and  x„;  the  time  dependency  is  indicated.  The 
values  of  u  computed  by  system.  (2.l4)  yield  a  sequence  u(t.), 
t„  <_  t.  £  t-,  5  which  allows  the  integrals  of  the  flux  over  the 
time  Interval  t-,-t„  to  be  approximated  more  accurately.   The 
sums  U-,  ,  U^   are  seen  to  be  telescoping  sums  in  f  cancelling  in 
pairs  over  all  net  points  between  x,  and  Xp. 

We  show  in  the  next  section  the  stability  properties  of 
system  (2.14);  indeed,  (2.l4)  is  unconditionally  unstable 
(after  all  this  effort!).   If  the  right  hand  side  of  (2.l4c) 
is  denoted  by  R  ,  then  a  stable  scheme  is  obtained  by  subtracting 
an  undivided  difference  quotient  of  fourth  order  from  R^ ,  i.e. 


(2.15) 


n+1    „n 

u.    =  R   - 


0)  ^4  n 
6  u. 


2? 


CO  >  0 


The  net  point  cluster  of  the  difference  scheme  (2.14)  is  shown 

below 

t 
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The  points  designated  °    are  first  order  accurate  data, 

•-  -point .3  are  second  order  accurate  data,  and  the  x_polnt  Is 

the  third  order  accurate  solution.   The  schematic  Indicates 

the  operator  defined  by  (2.l4a)  Is  applied  five  times; 

(2.l4b)  Is  applied  three  times,  and  (2.1^c)  and  (2.15)  are 

applied  once.   Another:"  scheme  Is  obtained  for  x,  =  2/3 

(ttpi-,  =  a^-.    =    1/3);  for  this  choice  of  t-,  the  difference  scheme 


Is  given  by 


(2.16a) 


(1) 


,n. 


^1+1/2  =  2  ^^1+1  +  ^l)  +  3  ^^^^1+1  -  ^1^^ 


(2) 


(1) 


,(1) 


(2  16b)   u^^^  .  u"^  +  -  {  ^  r(f^^^    -  f^^^   )  +  -(f""   -f"   )!} 
k^.iDo;   u^      ^1+3^2  L^  1+1/2     1-1/2^  ^  2^  1+1   i-l^^^ 


and  equations  (2.l4c)  and  (2.15). 

The  net  point  cluster  of  this  difference  scheme  Is  also 
shown  In  the  accompanying  figures. 


1-2 


1  +  2 


1-1     1     1+1 
FIRST  STEP 
(Five  applications  of  (2.l6a)) 


1-2    1-1     1     1+1    1+2 
SECOND  and  THIRD  STEPS 
(Three  applications  of  (2.l6b) 
and  one  application  of  (2.1^^0) 
and  (2.15)) 


Figure  3 
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( 2 ) 
Since  (2.l4c)  requires  second  order  data  u.   ,  as  well  as 

initial  data,  u.    for  evaluation,  in  principle  any  second  order 

difference  operator  can  be  used  to  generate  this  data;  in 

(2) 
particular  u.    could  be  obtained  from 

2 

u^^-*  =  u"^  -  -(f'^   -f*^  )+    ^^    (A      ff    -f  )-    A      (f  -f    )} 
^i      ^i    3^  i+1   1-1^    9  ^^i+1/2^  1+1   1^   ^i-1/2^  1   1-1^^ 


which  is  the  Lax-Wendroff  method. 

In  addition,  a  simple  method  can  be  used  to  obtain  first 
order  accurate  data  u.    other  than  (2.l4a)  or  (2.l6a).   The 
only  requirement  is  that  the  overall  scheme  be  stable.   In 
choosing  a^^   =    T-,   we  have  considered  only  values  which  result 
in  positive  weights  in  (2.5)  and  a  value  for  the  time  step  of 
(2.4)  which  is  less  than  or  equal  to  the  time  step  of  (2.5). 
The  permissible  range  of  a-,  ^  required  to  satisfy  these 
conditions  is  1/3  £  ct-i  n  f.  2/3. 

3 .   Stability  of  the  One  Dimensional  Difference  Operator 

Let  M  be  the  amplification  matrix  obtained  by  first  letting 
f(u)  =  Au  and  then  substituting  (2.lila)  and  (2.l4b)  into  (2.l4c) 
subject  to  the  viscosity  expression  (2.15).   The  amplification 
matrix  of  this  combined  system,  obtained  by  substituting 
Uq  exp  (iKx)  for  u"?  is 
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(3.1) 


2  2 
M(^,X,aj)  =  I  -  ^-|-  sln^?  -  ^  (1-cos  C)^ 


+  lAA  sin  ?(1  +  I  (1-cos  C)(1-X^A^)) 


where  ^  =  k  Ax  and  0  <_  C  1  f^  • 

Call  m  the  eigenvalues  of  M  (see  Figure  ^)  and  a  =  Xp, 
p  the  eigenvalues  of  A.   Construct  the  function  g(^,a,oj)  from 
the  real,  R,  and  imaginary,  1,  parts  of  (3.1)  by 


(3.2) 


2     2 
=  R   +  1   -  1 


=  Im^l  -  1 


Since  it  takes  so  much  effort  to  compute  g  we  shall  state  its 
value  here . 

g(C,a,aj)  =  ip  sm  5  +  j" (2  "*■  ^ o    )sin  ^(1-cos  ?) 


(3.2a) 


2 

+  ^(1-cos  ^)  -  ^(1-cos  ^)  +  T-  a  (1-a  )sin  ?(l-cos  5) 


Now  |m  I  <  1  if  and  only  if  g  <  0.   Allow  5  =  tt  and  observe  that 

(3.3)  g(TT,a,(jo)  =  ^  w(w-3)  <  0 

if  and  only  if  0  1  ^  £  3 . 

For  small  values  of  C,  m  can  be  written  as 


(3.^)      m(c,a,a3)  =  e^°^  +  (4a^-  a^-  o))  ^  +  0(5^) 
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Geometrically  one  sees  that  g  will  exceed  one,  in  the  complex 
plane,  unless 

(3.5)  0)  >  4a^  -  a   . 

That  is,  the  operator  (2.1^1)  will  not  be  stable.   Hence 

2   4 
combining  (3.3)  and  (3.5)  we  conclude  that  Ho   -o     ±  ^  ±  ^ 

2  2 

which  in  turn  implies  that  0<a   £l,   3£a   £.^.   To  show 

2 
that  a   takes  on  allowable  values  only  in  the  interval 

prescribed  by  the  Courant-Friedrichs-Lewy  condition  [4],  it  is 

necessary  to  show  that  for  any  value  of  o  such  that  3   ±  o      ±  H, 

there  exists  a   i   =    E,      such  that  g(5„,a,w)  >_  0  for  o)  given  by 

(3.5).   In  equation  (3.2a)  set  E,   =    tt/2: 


(^  ,a,a))  =  ^  (oj^+a)(6a^-12)  +  ( 4o^-23a^  +  28a^ ) ) 


2  '"''^^  -  JB 

>  ^  ((^a^-o^)^+  (4a^-a^(6a^-12) 

+  (ila^-23a^+28a^)) 

2 
=  £^  (a2-l)(a^-4)(a^-5)   1  0   for  3  <  a^  <  4 

2 
We  see  g  <  0  only  if  0  <_  a   £  1. 

,  2   4 
It  is  clear  that  if  oj  =  4a  -a   and  if  0  <  a  <  1,  then  g  <  0. 

2   i|  2 

Setting  gx(C,a)  =  g{E,,o,Ho   -o    ),   and  noting  sin  C  =(l-cos  ?) 


•(1+cos  5),  we  compute 

2 


2 
(3.6)    g*(5,a)  =  |z:  (1-cos  ?)2(4-a^) (l-a^) •P^(cos  C) 
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where  the  quadratic  polynomial 

2 
Pp(x)  =  ax   +  bx  +  c 

a  =  -  (l-a^) 
(3.7)  P 

b  =   2(3-a^) 

c  =  -  (S-a^) 

It  can  be  shown  that  Pp(x)  <  0  If  x  <  1,  which  Implies 

Pp(cos  C)  <  0  If  0  <  C  1  ■'T .   Clearly,  g^(^,a)  <  0  since 

2 
0  <  a   <  1.   One  also  observes  that  g^COjCT)  =  O*   At  this 

point  we  have  shown  that  the  right  hand  side  of  (3.2)  is 

negative  definite  for  tt  >_  C  >  0  and 

1)     0  <  a  <  1 
(3.8) 

ii)     00  =  4a   -  a 

We  exclude  the  case  where  a  =  1,  since  if  that  occurs,  co  =  3j 
and  g(^,l,3)  =  1,  which  leads  to  |m  |  =  1,  i.e.  the  associated 
difference  operator  is  not  dlsslpative. 

We  look  at  the  quantity  |m|.   It  can  be  bounded  (using  (3.^)) 
for  small  E,: 

(3.9)  |ml   <   1  +  (lla^-oa-aS  ^ipr  +    ^i^'^)    . 

One  can  show  that  for  any  a  e  (0,1)  there  exists  an  e  >  0  such 

,     2      H 
that  if  CO  =  4a  -a  +e   and  0  <  ^  <  it  ,  g(5,a,oo)  <  0,  which  is 

equivalent  to  |m|  <  1.   There  exists  a  pair  6-,,n  ,  each  greater 

than  zero  J  such  that  0  <_  C  £  n   implies  |ml  <_  1-6-,  C  •   Since 
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[fijTr]  is  compact  and  since  |  m  |  <  1  on  [n,TT],  there  exists  a 
positive  62  such  that  if  n  ;£  ?  £  ir,  then  |m|  <_  l-6pC  •   Let 
6  =  minimum  (6-,,6„);   we  see  that  the  difference  scheme  associated 
with  the  amplification  matrix  (3.1)  is  dlsslpatlve  for  0    <_  E,    <_  t\ 
in  the  sense  of  Kreiss  since 

(3.10)  |ml  <  1  -  6  5^  ,  6  >  0  . 

Since  the  accuracy  is  of  order  3,  (3.1)  is  stable  [5]. 

4 .     Two  Dimensional  Methods 

For  two  dimensional  hydrodynamic  flows  in  which  x  and  y  are 
the  cartesian  coordinates,  the  equations  of  motion  can  be 
written  in  conservation-law  form 

(^.1)  -t  =  ^x  -^  Sy 

where  g  is  the  vector  representing  the  flux  of  the  mass,  momentum 
and  energy  per  unit  volume  in  the  y  direction.   We  carry  out 
differentiation  of  (4.1)  using  the  chain  rule  to  obtain 

(4.2)  u,  =  A(u)u   +  B(u)u   . 

■c         x         y 

In  general  the  miatrices  A  and  B  do  not  commute  and  are  not  normal. 

If  one  considers  the  class  of  linear  problems  where  A„  =  A(u„) 
and  B   =  B(Uj.),  u„  the  state  about  which  the  motion  Is  linearized, 
then  (4.2)  may  be  Integrated  to  yield 

u(t  +  At)  =  P  u(t)  , 
(4.3) 

P  =  exp  [(Aq  |3^+  Bq  |^)At]  . 
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Equation  (4.3)  is  also  valid  even  if  A^  and  B^  vary;  however, 
the  variation  must  be  independent  of  time.   With  the  obvious 
change  in  notation  P  may  be  written  as 

(4.4)  P  =  e^"^^  . 

This  operator  is  called  the  exact  solution  operator  of 

equation  (4.2).   The  Fourier  transform  of  P,   P  =  exp  (i(AC+Bn)) 

is  called  the  symbol  of  the  operator   P    (see  ref.  [5]). 

By  multiplying  the  initial  data     u(x,y,0)  =   u(0),  r  successive 

times  using  the  operator  P,  we  can  map  u(0)  into  u(T),  T  =  r  At . 

In  forming  difference  approximations  to  (4,1)  or  equivalently 
in  approximating  the  operator  P,  the  question  of  stability  arises. 
The  analysis  of  stability  of  the  difference  operator  becomes 
difficult,   especially  as  the  order  of  accuracy  (and  corresponding 
complexity)  of  the  difference  scheme  increases.   Indeed,  in 
Strang's  paper  [6]  on  the  construction  of  accurate  difference 
methods,  he  is  motivated  in  the  construction  of  difference  methods 
by  approximating  P  to  the  desired  degree  of  accuracy. 

It  is  well  known,  for  instance,  that  a  first  order  approxima- 
tion to  the  matrix  P  can  be  written  as 

(4.5)  P  =  e^e^  +  O(At^) 

since  A  and  B  are  of  order  At.   The  error  in  (4.5)  goes  to  zero 

A       B 
when  A  and  B  are  scalars.   The  operators  e   and  e    can  be 

thought  of  as  exact  solution  operators  to  the  one  dimensional 

differential  equation  of  the  form  (2.2a)   defined  separately 
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for  the  X  and  y  directions.   Let  L(A)  and  L(B)  be  difference 

A        R 

approximations  to  the  operators  e   and  e   respectively.   If 
L(A)  is  the  symbol  of  the  x-difference  operator  and 

(4.6)  |L(A)  -  e^^^l   -  0(AtP"^^) 

then  we  conclude  that  the  difference  operator  is  accurate  to 
order  p . 

Strang  has  shown  [6,7,81  that  if  one  considers  an  operator 
L-,  (A,B)  formed  from  the  product  of  the  one  dimensional  operators 

L^(A,B)  =  L(A)  L(B) 
then 

(4.7)  ||  (L^(A,B)  +  L^(B,A))  -  P|  =  O(^^n^) 

where  P  has  been  defined  previously.   Strang  has  also  noted 
recently  [8]  that  it  is  possible  to  satisfy  (4.7)  with  the  product 

(4.8)  L(A/2)  L(B)  L(A/2)  , 

replacing  the  sum  in  (4.7);  hence   (4.8)  provides  the  structure 
for  another  difference  scheme  of  second  order  accuracy. 

The  stability  of  (4.8)  follows  immediately  from  the  stability 
of  each  one  dimensional  operator  L(A)  and  L(B).  For  (4.8)  to  be 
second  order  accurate  and  stable  each  one  dimensional  operator, 

given  by 

2 

(4.9)  L(A/2)  =  I  -  I  A6  +  ^  A^6^  , 

which  is  the  Lax-Wendroff  operator,  need  be  stable;  this 
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requirement  is  fulfilled  if  the  eigenvalues  of  A  and  B,  y(A) 
and  ij(B)   satisfy 

X  y(A)  <  1 
(4.10) 

X  y(B)  <  1  . 

Equation  (4.8)  can  be  used  with  a  second  order  two  step  procedure 
rather  than  (4.9),  i.e.  equations  (2.l4a)  and  (2.l4b)  or  system 
(2.16)  (with  the  appropriate  time  step).   The  advantage  is  the 
elimination  of  the  evaluation  of  the  matrix  A  in  the  difference 
scheme . 

Gourlay  and  Morris  [9]   have  performed  some  computations 
with  such  schemes.   They  have  adopted  the  operator  in  (4.7) 
for  practical  computations  by  also  using  two  step  versions  of 
L(A)  and  L(B) . 

We  wish  to  look  for  difference  schemes  of  the  form  given  by 
(4.8)  which  are  of  uniform  third  order  accuracy.   The  structure 
of  the  difference  operator  will  then  be  based  on  the  third  order 
one  dimensional  operators  discussed  in  Section  2. 

We  have  considered  generalizations  of  the  operators  given 
in  (4.7)  and  (4.8)  of  the  form 

a.  .A   6.  .B 

(4.11)  S   =    I    c.    ri  e   ^^      e  ^J   ,        I  c   =  1 

J   -^   i  J   ^ 

and      y  a. .  =  y  6. .  =  1  with  each  a, 3  >  0. 

Clearly  if  one  chosses  c,  =  1  with  ct-, -,  =  3-,-,  =  gpi  ~  "^qi  ~  ■'-/^ 
and  api  =  0  =  3o-|   then  (4.11)  becomes 

(4.12)  S^(A,B)  =  e^/2  e^  e^^^  . 
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With  the  constants  c,  =  1/2  and  a-,-,  =  3-,-,  =  1,  and  Cp  =  1/2 
and  a,  p  =  0  =  Ppp,  g-,p  =  1  =  °'22'  (^-H)  becomes 

(4.13)  ^2  =  I  (e^e^  +  e^e^)  . 

The  operators  S-,  and  Sp  are  the  operators  Strang  has  investigated 
It  is  interesting  to  consider  the  operator  formed  from  linear 
combinations  of  (4.12)  and  (4.13),  i.e. 

.    S, (A,B)+S, (B,A)     . 

(4.14)  S^  =  ^(^ 2"^ )  -  ^  Sp  . 

If  the  one  dimensional  differential  operators  in  (4.l4)  are 
replaced  by  corresponding  one-dimensional  difference  operators 
defined  by 

2         3 

(4.15)  L(A)  =  I  +  AAy6  +  |-  A^6^  +  ^  A^fi^yfi   , 

then  it  may  be  verified  by  direct  computation  that  the  resulting 
difference  approximation  L^  to  the  differential  operator  S^ 
satisfies 

2  3 

(4.16)  L^  =  I  +  A(A+B)  +  ^  (A+B)^  +  j|-  (A+B)^  . 

Here,  for  simplicity  we  have  used  the  abbreviations 

A  =  A(u)6  ,   B  =  B(u)6    for  the  centered  difference  operators 

^         y 

and  X    =    At/Ax,   A  =  Ax  =  Ay .   This  is  precisely  the  expansion 
for  (4.4)  up  to  cubic  terms.   This  operator  L-,  was  first  found 
by  J.  Dunn. 

The  procedure  used  to  derive  the  third  order  approximation 
to  (4.4)  which  is  in  some  sense  computationally  optimal  follows. 
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The  operator  defined  by  (4.l4)  is  complicated  and  inefficient 
as  it  requires  ten  sweeps  through  the  mesh  —  five  in  the 
x-directlon  and  five  in  the  y-direction  —  to  advance  the 
solution  one  time  step.   It  is  clear  that  more  compact  forms 
resulting  in  economical  algorithms  suggested  by  (^.11)  are 
desirable.   Consider  the  differential  operator 

M,  -l'7^  A  B  ^     aA  BB  (l-a)A  (1-3)B 

(4.17)         c-,e  e   +  Cpe   e   e      e 

If  the  constants  are  chosen  correctly  (^.17)  can  be  made 
to  differ  from  (^.4)  by  terms  of  0(At  ).   To  do  this  first 
expand  each  of  the  exponential  forms  up  to  terms  involving 
cubic  powers  of  the  matrices  A  and  B.   This  allows  the  evaluation 
of  the  term  corresponding  to  the  coefficient  c^ 

I  +  (A+B)  +  I  (A^+B^)  +  [a+(l-a)(l-B)AB 

+[6Cl-a)]BA  +  ^(A^+b3)+  |[a^+2a(l-a)(l-B)+(l-a)^(l-6)]A^B 
(4.18) 

+  a6(l-a)ABA  +  |  6(l-a)^BA^  +  |  6^(l-a)B^A 

+  3(l-a)(l-B)BAB  +  |[aB^+  2a6(l-6)  +  (1-6)^]AB^  , 

and  the  evaluation  of  the  term  corresponding  to  the  coefficient  c, 

(4.19)    I  +  (A+B)  +  |(a2+b2)  +  AB  +  |(a3+b3)  +  ^(A^B  +  AB^)  . 

For  simplicity  we  just  equate  the  coefficients  of  the  matrices 
BA,  ABA  and  BAB  in  (4.l8)  to  their  proper  values: 
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and 


c^BCl-a)   =   2   ' 
C2a3(l-a)  =   ^   5 


C26(l-a)(l-B)  =  I 


These  equations  yield  the  values  a  =  1/3,   3  =  2/3,   c^  =  9/8 
and  C-,  =  -1/8.   We  again  observe  the  appearance  of  the  nonpositive 
weight  in  our  difference  scheme.   The  difference  equation  becomes 

(it. 20)  S^  =  I  L(A/3)L(2B/3)L(2A/3)L(B/3)  -  ^  L(A)L(B)  . 

Each  one  dimensional  difference  operator  in  (4.20)  is  defined 
by  (4.15).   The  proof  of   stability  of  (4.20)  (and  that  of  (4.l4)) 
does  not  follow  from  the  fact  that  the  norm  of  each  one 
dimensional  operator   |l|  <   1.   If  each  coefficient  c.  in  (4.17) 
were  greater  than  zero  then  S-.  would  be  a  convex  operator  and 
one  could  conclude  in  that  case  that  S^  had  norm  less  than  one. 
We  defer  this  question  until  later. 

It  appears  that  (4.20)  is  most  efficient  in  the  sense  that 
the  number  of  one  dimensional  sweeps  is  a  minimum  for  a  third 

order  operator.   One  needs  at  least  six  applications  of  the 

A       B 
exponential  operators  e   and  e   to  match  the  noncommutative 

terms  that  result  from  the  third  order  term  in  the  Taylor 

expansion  for  e    ;  i.e.  z-(A+B)  .   The  proof  of  this  statement 

involves  consideration  of  linear  combinations  of  products  of 

A      B 
e   and  e   taken  two  at  a  time  and  three  at  a  time.   All  such 

combinations  fail  to  yield  simultaneously  the  matrix  operators 

ABA  and  BAB.   Next  consider  product  combinations  of  the  one 
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dimensional  operators  taken  four  at  a  time: 

aA   BB  ^(l-a)A  ^(1-6)B 
e    e   e       e 

Expanding  and  considering  the  requirement  that  second  order 

accuracy  Implies  B(l-a)  =  1/2,  we  find  that  for  third  order 

accuracy  a  =  1/3,  3  =  2/3.   Hence  there  Is  a  contradiction. 

Finally  operators  formed  from  products  taken  five  at  a  time 

are  of  the  general  form 

a-,A   B-|B   ttpA   BpB   a^A 
(4.21)  e^   e^  e    ^      e    e-^   . 

For  (4.21)  to  be  third  order  accurate  a^  must  satisfy 

2 
12a-,  -  6a-,  +1  =  0.   This  polynomial  however  has  only  complex 

roots . 

Hence  a  third  order  splitting  method  of  the  form  (4.11) 
must  have  at  least  six  terms.   We  state  that  the  linear  combina- 
tion of  the  form  (4.12),  I.e. 

aA  3  (l-a)A  ,     3B  A^(1-3)B 
C-,  e   e  e       +  c^e   e  e 

cannot  differ  from  e     by  terms  of  0{^t    ).   Satisfying 
consistency  requires  that  a  and  3  must  satisfy 

|(a-3) (a+3-1)  =  0  . 

If  a  =  3  we  can  show  that  a  must  satisfy  a  quadratic  In  a  with 
complex  roots.   If  a+3  =  1,  a  satisfies  (a-l/2)(a-I)  =  0. 
Now  a  7^  I  so  a  =  1/2  and  therefore  3  =  1/2  which  Is  a  contradiction, 
If  one  were  to  only  consider  operators  with  c.  >  0,  then 
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third  order  accuracy  could  still  be  obtained  but  with  a  relaxa- 
tion of  the  condition  that  all  a.,  B.  >  0.   Although  the  analysis 
of  stability  would  be  trivial,  one  would  have  to  accept  multistep 
difference  methods  with  operators  having  a  negative  time  step. 
For  flows  which  contain  shocks  or  other  irreversible  phenomena 
the  problem  is  not  well  posed.   If  the  flow  is  smooth  and 
thermodynamically  reversible  there  may  be  no  drawback  to  such 
methods.   We  indicate  in  section  7  some  results  using  (4.20). 

5 .   Asymptotic  Operators 

It  is  possible  to  generate  a  positive  difference  operator  but 
only  asymptotically.   Consider  the  differential  operator 

S(A,B;N)  E  e^/N  ^2B/N  ^A/N 


and  its  conjugate  S  (B,A;N).   Then 


C5.1) 


N/4 


S„  =  (S(A,B;N)   S(B,A;N))     ,    N  -  4,8,... 


N 


is  called  an  asymptotic  third  order  difference  operator.   S 


N 


wo 


uld  be  an  exact  third  order  operator  if  the  coefficient   6-, 
of  the  terms  (A^B,B^A),  (BA^,AB^)  and  the  coefficient  62   of 
the  terms  (ABA,  BAB)    satisfied  6-.    =    6^   -    1/6.   Instead  these 
coefficients  are  functions  of  N.   We  have  computed  bounds  on 
the  coefficients  and  show  them  below  for  several  values  of  N: 


N 


12 


I  61-1/6  I  < 


.0053 
.0013 
.0006 


62-1/6  I  < 


ni05 
0026 

0012 
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It  appears,  that   using  (4.4), 

|llm  S   -  P|  =  0(AtS  . 

Since  |Sj^(A,B)|  <  1  and  |Sj^(B,A)|  <  1,  then  |  S^  |  <  1   which 
shows  the  stability  of  (5.1).    The  operator  defined  by  (5.1) 
achieves  its  accuracy  by  using  finer  and  finer   time  steps, 
At/N,  as  N  ->-  ~;  indeed  it  is  the  form  assumed  for  S„  that  gives 

the  rapid  convergence  of  coefficients  &-.  ,    6p.   The  operator 

/  A/N  B/NnN/2   m    o  )i        -11       ->  ■  ■  ^-   u-  u 

(e  '  e    )    ,  N  =  2,4,...   will  also  give  asympcotic  high 

order  accuracy  (greater  than  first  order)  but  requires  many 

evaluations  (large  N)  per  time  step.   In  comparison,  (4.22)  may 

be  satisfactory  for  N  =  4. 

6 .   Stability   of  Two  Dimensional  Operators 

Except  for  the  brief  discussion  on  the  stability  of  asymptotic 
operators,   we  have  not  found  a  satisfactory  method  for  the 
analysis  of  the  stability  of  the  operators  given  by  (4.11).   Our 
only  recourse  is  to  carry  out  a  numerical  analysis  of  the  eigen- 
values of  the  amplification  matrix  using  the  digital  computer. 
We  have  completed  a  calculation  in  which  the  independent  variables 

are  the  wave  numbers  (^,ri)  in  (x,y)  space.   We  took  the  dissipa- 

,  2   4 
tion  coefficient   to   to  be    oj  ==  4a  -a  +e  with  -0.2  ±  e   ±   0.2  in 

steps  of  0.1  and  with  0  <_  a  <_  1  also  in  steps  of  0.1. 

Prom  this  parametric  study,  it  appears  that  the  spectral 

radius  of  S^iE,,n),    the  transform  of  (4.20),  satisfies 


-26- 


(6.1)  lyCs^CC^n)) I  <  1 
and 

(6.2)  |y(S2(?,n))  I  <  1  -  6|0l^ 

2  2    1/2 
where  6  =  (5  +n  )     is  the  L„  norm  in  wave  number  space  If 


1)       0<a<l,        £>0 
(6.3) 

ii)        a   ==  0   ,        e  >  0  . 

However  if  0  <  a  <  1  and  e  <  0   the  spectral  radius  exceeds  one 
Indeed  if  e  =0.1   we  can  choose^  in  (6.2),  6  =  10    uniformly 
Independent  of  a.   To  achieve  this  define 

e 

then  pick 

6  =   inf   6 '  . 
(C,n) 

We  have  found  that  6'  is  smallest  when  (^,n)  is  near  (tt,tt). 
For  5  =  n  we  indicate  the  behavior  of  |  y  (  E;  ,  C )  |  by  the 
following  table. 
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r 

( 

yCS^CCO)    -1) 

K     =     -.1 

K    =     0 

K     =     +.1 

0.0 

0.0 

0.0 

0.0 

0.01 

2.9x10"^^ 

-1.6x10"-^^ 

-6.0x10"^^ 

0.10 

2.9x10""^ 

-1.6x10"'^ 

-6.0x10"''' 

0.5 

1.7x10 

-9.7x10"^ 

h 
-3.7x10 

1.0 

2.6x10"^ 

-1.5x10"^ 

-6.5x10"^ 

2.0 

3.2x10"^ 

-1.9x10"^ 

-6.7x10"^ 

3.0 

1.0x10"^ 

-6.2x10"^ 

-2.1x10"^ 

We  have  defined  the  distance  from  the  origin  r  =  /2  ^  in 

2   4 
wave  number  space  and  k  =  oo-(4a  -a  )  .   In  the  next  section  we 

present  further  evidence  as  to  the  usefulness  of  these  third 

order  operators . 

7 .   Results 

We  describe  some  numerical  experiments  carried  out  with  the 
scheme  (2.14)  for  the  Riemiann  problem  in  one  dimension  and 
with  (4.20)  for  a  two  dimensional  scalar  problem  invented  by 
Crowley  [10]. 

Figures  5  and  6   show  the  results  of  two  calculations 
using  system  (2.l4)  and  (2.15)  to  obtain  approximate  solutions 
to  (2.2).   Both  calculations  start  with  the  same  initial  data, 
i.e.  two  constant  states  separated  by  a  discontinuity: 
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p"^(x) 
u(x) 
p(x) 


2 

0 

0 

0 

0 

.571 

'X    >    0 


p"^(x) 
u(x) 
p(x) 


2.245 
0.698 
3.528 


X  <  0 


In  Figure  5   a  =  0.9  with  to  =  0.75.   The  instability  is  clearly 
shown.   Figure  6  shows  the  solution  at  approxinately  the  same 
time;  here  a  =  0.9  with  to  =  2.5. 

We  have  tested  the  third  order  m.ethod  in  two  dimensions  on 
the  following  scalar  problem.   The  differential  equation 


(7.1) 


r,  +  ur   +  vr   =0 
t      X     y 


describes  the  motion  of  the  function  r(x,y,t)  in  the  x-y  plane 
if  the  velocity  components  u  and  v  are  specified.   We  take  them 
to  be 
(7.2) 


u 

V 


-y 

X 


which  means  that  the  velocity  vector  depends  only  on  the  radius, 
i.e.  V  =  fi  r   defines  solid  body  rotation  (in  our  problem 
centered  at  (x^jy^)  =  (30,30)).   If  the  components  of  (7.2)  are 
differentiated  with  respect  to  x  and  y  respectively  we  see  that 
one  may  write  (7.1)  in  conservation  form 


(7.3) 


r,  +  (ur)^  +  (vr)   =  0 
L-        A        y 


since  the  velocity  field  is  divergence  free. 

The  distribution   r(x,y,t)   is  prescribed  at  t  =  0  to  be 

a  right  circular  cone  in  (r,x,y)  space  centered  at  (37,37) 

with  base  radius  of  five.  Ax  =  Ay  =  1 .   Equation  (7.3)  subject 
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to  (7.2)  states  that  the  total  variation  of  r(t),  Dr(t)/Dt, 
along  circles  with  radius  centered  at  (x„,y„),  vanishes,  i.e. 
r  is  just  uniformly  rotated  with  period  t  =  2TT/fi.  In  our 
computation  the  mesh  size  is   60  x  60   while  the  cone  height 
is  one. 

The  table  shown  below  is  a  summary  of  the  computations 
performed  for  this  "cone"  problem.   In  problem  1,  the  first 
order  scheme  is  defined  by  the  operator  (4.5)  while  for  problems 
2,  3  and  4,  the  second  order  scheme  is  defined  by  the  operator 
(4.8).   The  third  order  scheme  is  given  by  (4.20)  with  L  defined 
by  system  (2.l4)  and  (2.15).   The  value  of  to  in  problem  5  did 
not  satisfy  the  stability  condition  (3.5);  it  was  kept  constant. 
For  problems  6,    7  and  8,  the  local  value  of  oj  satisfies 
to  =  4a   -  a   +  £  with  e    =  .01.   The  value  of  a  =  At/A  |u|  with 
Iu|  =  (u  +v  )    .   The  components  of  drift  of  the  vertex  of  the 
cone  in  the  x(y)  direction  equals  the  x(y)  position  of  the  vertex 
computed  by  the  difference  method   -  the  x(y)  position  of  the 
vertex  given  by  the  exact  solution. 

We  see  how  poorly  first  order  methods  compare  with   second 
or  third  order  methods  in  the  amplitude  and  phase  of  the 
solution.   The  most  striking  difference  between  second  and 
third  order  accuracy  is  in  the  computation  of  the  phase  of  the 
solution.   The  position  of  the  vertex  is  within  one  half  mesh 
width  in  both  the  x  and  y  direction  for  the  third  order 
calculation   but  is  two  to  three  mesh  widths  from  the  exact 
position  in  both  the  x  and  y  direction   for  the  second  order 
calculation.   For  both  second  and  third  order  schemes,  increasing 
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Table.   Summary  of  Computations  for  Cone  Problem. 


r-   '■- -"I 

Computed 

Vertex 

Drift 

No.  of 

Rotations 

Vertex 

- 

Integra 

Problem 

Method 

Traversed 

Amplitude 

x-direction 

y-directlon 

-tion 
Cycles 

1 

first  order 

1/^ 

.07856 

6.28699 

2.07856 

150 

2 

second  ordej 

1 

.98935 

1.65263 

-2.51629 

600 

3 

second  ordei 

1 

•98363 

1.55743 

-2.29384 

300 

%ax=  1/3 

2 

.82304 

2.37585 

-3.70616 

600 

H 

second  orde: 

2 

.79365 

2.25867 

-3.41507 

240 

a   =  5/6 
max 

1 

5 

third  order 
w=const  .  = . O; 

1 

1.15205 

.33053 

f 
1 

-  .29569 

600 

6 

.third  order 

'a        =  1/6 
max 

1 

1.03803 

.26070 

-  .34469 

600 

7 

^third  order 

1 

.99707 

-  .24833 

-  .36810 

300 

%ax=  1/3 

2 

.89400 

-  .89400 

-  .45386 

600 

8 

third  order 

a   =  ^/7 
max 

O 

£_ 

.81353 

-  .33365 

-  .45469 

350 

9 

third  order 

a        =  ^/6 
max   ^ 

^  2 

unstable 

10 

third  order 

5/12 

unstabJ.e 

11 

third  order 

1/2 

goes 

150 

i)=oonst  .  =  .  01 

1 

unstable 

300 

Wx  =  1/3 
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the  time  step  At  with  fixed  space  step  increases  a  which 
results  in  greater  dissipation  in  the  third  order  difference 
scheme.   Increasing  At  also  increases  the  artificial  viscosity 
in  the  second  crder  method  [see  [5]]  and  therefore  the  greater 
smoothing  reduces  the  maximum  amplitude  of  r(x,y,t). 

Problems  9  and  10  went  unstable  for  the  values  of  a 
indicated.   Hence,  one  obtains  an  approximate  upper  bound  for  o, 
which  gives  an  approximate  upper  bound  for  an  allowable  time  step. 

The  remaining  figures  are  labeled  as  to  problem  number,  which 
corresponds  to  the  problems  given  in  the  table  on  the  preceding 
page.   The  figures  show  the  overall  behavior  of  the  various 
methods  and  give  means  for  a  quick  comparison  between  the  methods. 
The  contour  lines,  at  each  instant  of  time,  define  values 
r(x,y)  =  constant,  the  values  of  which  lie  between  0.05  and  0.95. 
For  clarity  the  snapshot  of  the  solution  at  the  latest  time  has 
been  shifted  by  an  amount  D  along  a  line  connecting  the  center  of 
rotation    and  the  vertex  of  the  cone. 

The  scheme  (^.20)  required  approximately  4  seconds  per 
sweep  while  the  second  order  method  (^.8)  (alternate  sweeps  were 
computed  using  first  L(A/2)  L(B)  L(A/2)  then  L(B/2)  L(A)  L(B/2), 
etc.,  rather  than  L(A/2)  L(B)  L(A)...  L(A)  L(B)  L(A/2))   required 
approximately  1  1/3  seconds.   By  comparing  the  numerical  results 
in  the  above  table,  it  appears  that  the  mesh  ratio  for  third  order 
methods  can  be  increased  by  a  factor  of  three  over  the  second 
order  method.   Comparable  errors  in  the  amplitude  of  the  solution 
are  obtained  with  the  two  methods  but  a  clear  superiority  in  the 
phase  of  the  solution  is  achieved  with  (4,20). 
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Our  tentative  conclusion,  subject  to  additional  numerical 
tests  is  that  (^.20),  using  a  more  coarse  mesh,  may  be  as 
economical  as  a  second  order  calculation  on  a  fine  mesh  while 
still  giving  superior  numerical  results. 
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Figure  4 

Absolute  value  of  the  eigenvalues  of  Equation  (3.1) 
showing  dependence  on  u.  The  values  of  o  are  within  ,05 
of  the  maximum  allowable  for  each  co. 
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Figure  5 

Pressure  and  (iensity  profiles  for  the  Rlemann 
problem  after  l88  At  (t  =30.138)  using  o  and  w 
not  satisfying  Equatic^n  (3-5). 
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Flgvire  6 

Pressure  and  density  profiles  for  the  Riemann 
problem  after  I63  At  (t=30.065)  using  a    and  w 
satisfying  Equation  (3.5).   The  rarefaction  wave 
propagates  to  the  left;  the  contact  discontinuity 
is  located  at  x  =  50;  the  shock  propagates  to 
the  right  with  an  error  less  than  1%    of  the 
theoretical  shock  speed. 
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Figure  7 

Problem  2  -   Second  order  method  with  a     =  1/6  and 

max 

computed  vertex  amplitude  equal  to  O.989; 
the  exact  value  is  1.0. 
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Figure  8 

Problem  3  -   Same  Initial  data  and  method  as  in  Figure  7 

but  with  a     =  1/3;  after  300  cycles  computed 

max       ^  J  i- 

vertex  amplitude  equals  O.983;  after  6OO  cycles 
amplitude  equals  0.823. 
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Figure  9 

Problem  6  -   Same  Initial  data  as  In  Figure  7;  third  order 

method  with  o^^^   =    1/6;   co  Is  variable  and  Is 
max 

computed  from  Equation  (3.5).   The  amplitude 
Is  1.03  after  600  cycles. 
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Figure  10 


Problem  7  -   Same  as  Figure  9  but  with  o^^^   =  1/3 


The  amplitude  after  300  cycles  is  0.997  and 
after  600  cycles  it  is  0.894. 
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Figure  11 

Problem  11  -   Same  initial  data  as  in  Figure  7;  third  order 

method  with  co  =  0.01  and  a     =  1/3. 

max 

Calculation  does  not  satisfy  stability 
condition  (3.5).   Eddies  are  forming  while 
the  amplitude  increases  —  calculation 
eventually  goes  unstable. 


•50- 


6.5 


5.5 


4.5- 


3.5- 


2.5 


.5- 


0.5 


d 

1 


Iteration  number  300 


Iteration  number  150 


__l 1 I I I I I I I I ! I 

0.5  1.5  2.5  3.5  4.5  5.5  6.5 


-51- 


This  report  was  prepared  as  an  account  of 
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person  acting  on  behalf  of  the  Commission: 
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the  use  of  any  information,  apparatus, 
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As  used  in  the  above,  "person  acting  on  behalf 
of  the  Commission"  includes  any  employee  or 
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such  contractor,  to  the  extent  that  such  em- 
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mation pursuant  to  his  employment  or  contract 
with  the  Commission,  or  his  emplojrment  with 
such  contractor. 
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